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ABSTRACT 

In recent work we presented the first results of global general relativistic magnetohydrodynamic 
(GRMHD) simulations of tilted (or misaligned) accretion disks around rotating black holes. The sim- 
ulated tilted disks showed dramatic differences from comparable untilted disks, such as asymmetrical 
accretion onto the hole through opposing "plunging streams" and global precession of the disk powered 
by a torque provided by the black hole. However, those simulations used a traditional spherical-polar 
grid that was purposefully underresolved along the pole, which prevented us from assessing the be- 
havior of any jets that may have been associated with the tilted disks. To address this shortcoming 
we have added a block-structured "cubed-sphere" grid option to the Cosmos-f-|- GRMHD code, which 
will allow us to simultaneously resolve the disk and polar regions. Here we present our implementation 
of this grid and the results of a small suite of validation tests intended to demonstrate that the new 
grid performs as expected. The most important test in this work is a comparison of identical tilted 
disks, one evolved using our spherical-polar grid and the other with the cubed-sphere grid. We also 
demonstrate an interesting dependence of the early-time evolution of our disks on their orientation 
with respect to the grid alignment. This dependence arises from the differing treatment of current 
sheets within the disks, especially whether they are aligned with symmetry planes of the grid or not. 
Subject headings: accretion, accretion disks — black hole physics — methods: numerical — MHD — 
relativity 
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1. INTRODUCTION 

We have recently undertaken a series of numerical 
studies of titled accretion disks around rapidly rotating 
black hol es, first in the hydrodynamic (Fragile & Anni- 
nos|2005|) and then in the magnetohydrodynamic (MHD) 
( [ii^'ragile et al. 2007b I limits. All of these simulations have 
been fully general relativistic, using the Kerr-Schild met- 
ric to represent the spacetime of the black hole. 

Tilted accretion disks are of particular interest because 
they are subject to differential warping due to the Lense- 
Thirring precession of the rotating black hole. For very 
thin disks, close to the black hole the competition be- 
tween the differential twisting and "viscous" damping 
causes the angular momenta of the disk and hole to align. 
Further out in the disk, beyond some warp radius, the 
disk maintains its misaligned state. 

For moderately thin to thick disks, such as the ones 
we simulated previously, the situation is more complex 
and interesting. The primary difference is that warp- 
ing is transported via bending waves rather than dif- 
fusively, as for thin disks. One consequence of this is 
that the midplane of a thick disk does not tend to align 
with the symmetry plane of the black hole at small radii, 
as in the thin disk case. In fact, the relative tilt be- 
tween the black-hole and disk angular momenta can in- 
crease at small radii. Having the tilted disk penetrate 
very close to the black-hole has many interesting conse- 



quences. For instance, we found that accretion onto the 
hole occurs predominantly through two opposing "plung- 
ing streams" that start from high latitudes with respect 
to bot h the black-hole and disk midplanes (Fragile et al. 
2007a). There is also a strong epicyclic driving within 
the disk attributable to the gra vitomagnetic torque o f 
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the misaligned (tilted) black hole ( Fragile fc Blaes|2008 1 
The induced motion of the gas can be coherent over the 
scale of the entire disk. The gas also experiences peri- 
odic (twice per orbit) compressions. The compressions 
occur as the gas orbits past the line-of-nodes between 
the black-hole symmetry plane and disk midplane. Near 
the black hole these compressive motions can become 
supersonic and transform into a pair of quasi-stationary 
shocks. The shocks act to enhance angular momentum 
transport and dissipation near the hole, forcing some ma- 
terial to plunge toward the black hole from well outside 
the innermost stable circular orbit. Finally, because we 
are simulating disks with finite radial extents and fast 
sound-crossing times, the torque of the black hole causes 
the entire disk body to precess globally. 

The main shortcoming of our work so far comes from 
limitations imposed on us by our use of a spherical-polar 
grid. First, construction of a uniform spherical polar 
grid in three-dimensions results in very small zones sur- 
rounding the two poles, where all of the lines of longitude 
converge. These very small zones constrain the Courant- 
limited timestep to be exceedingly small, such that the 
required CPU cycle count becomes prohibitively large. 
To avoid this problem, researchers have eith er excised a 
small conical section around each pole (e.g. De Villiers 
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fc Hawleyll200 3j or u sed a lower grid resolution near the 



poles (Fragile et al. 2007b I. Although these techniques 



are reasonable when one is primarily interested in study- 
ing the equatorial region (where a disk may form) , these 
are not satisfactory when one is interested in what is 
happening in the polar regions (where jets may form). A 
second problem with the spherical-polar grid is that the 
poles themselves actually represent coordinate singular- 
ities, which present significant challenges for numerical 
advection and curvature coupling schemes (e.g. solving 
Riemann curvature source terms). 

For these reasons we have added t he cubed-sphere grid 
(Sadourny 1972 Ronchi et al.||1996" l as an option within 
our numerical code, Cosmos-|--|-. 'ihe advantage of this 
grid construction is that its topological properties more 
closely resemble a Cartesian coordinate system than a 
spherical-polar system. The cubed-sphere grid uses a 
more uniform zone spacing than spherical polar, so the 
timestep can remain reasonably large even in high res- 
olution simulations. Also important, the grid does not 
contain any coordinate singularities except at the ori- 
gin, which is not a concern for our intended use since 
we truncate the grid just inside the event horizon of the 
black hole. Ours is not the first apphcation of the cubed- 
sphere grid to problems in computational astrophysics; it 
has been used previously to study a ccretion onto rotating 
stars with inclined magnetic fields ( Koldoba et al.|[2002 



Romanova et al.|2003[) a nd a few pr oblems in stellar evo 



lution (Dearborn et al. 112005 20061. However this is the 



first application of this grid to the study of black-hole 
accretion disks and their attendant jets. 

The paper is organized as follows: Sj2] describes the 
cubed-sphere mesh in detail and our particular imple- 
mentation. In S|3]we discuss results of basic gradient tests 
on the cubed-sphere grid. In S|4]we compare two sets of 
numerical simulations of black-hole accretion disks. In 
the first set we compare simulations of disks accreting 
onto a Schwarzschild black hole. We compare differ- 
ent grids, resolutions, and orientations of the disk with 
respect to the grid. Because these simulations use a 
Schwarzschild black hole, the orientation should have no 
physical meaning. However, we show that there are, nev- 
ertheless, considerable differences in their evolution at 
early times. We present the case that the differences 
have to do with the differing treatments of the midplane 
current sheets in the disks, which forms from the differ- 
ential winding of our initial poloidal field loops. Finally 
we compare simulations of tilted accretion disks around 
Kerr black holes. We use a tilt of /?o = 15° and a spin of 
o/Mbh — t/BH/-^BH = O-^i ™ geometrized units where 
G = c = 1 , and Mbh and Jbh are the mass and angular 
momentum of the black hole, respectively. One of these 
simulations is run on a spherical-polar mesh, the others 
on the cubed-sphere grid. We demonstrate that, for the 
most part, the simulations agree very well. 

2. THE CUBED SPHERE 

The cubed-sphere grid gets its name from its construc- 
tion - it is actually composed of six "blocks" that are 
morphed into segments of a sphere. Each block is con- 
structed of segments of concentric radial shells. In the 
present work, these shells are spaced exponentially based 
upon their distance from the hole, similar to a logarith- 
mic radial coordinate. The other two coordinates are 



constructed such that, on any given block, the grid lines 
trace out "great circles" on each radial shell segment. It 
is as if there are two longitude coordinates, 0i and 4)2, 
on each block. The range of and 02 on each block is 
7r/2 so that the full Air steradian is covered by the six 
blocks. 

The difficulty with the cubed-sphere grid is that the 
"great circles" cannot be made continuous across all six 
blocks, and hence the block-structured nature of the 
mesh. Stated differently, the coordinates 0i and 4)2 can- 
not maintain a consistent orientation across all blocks. 
At each block boundary, the coordinate system has a 
discontinuous jump. Fortunately this can be handled 
with the proper application of boundary conditions and 
communication between blocks, as we shall describe. 

Another problem with the cubed-sphere grid is that 
the 4i and 4^2 coordinates are not orthogonal. There 
are techniques available to try to improve the orthog- 
onality of the cubed-sphere grid at the cost of reduc- 
ing its uniformity. However, such techniques have been 
shown not to perform significantly better th an the stan- 
dard cub ed-sphere grid implemented here (Putman & 
Lin 20071. Furthermore, such techniques are not nec- 



essary m our Cosmos-|— I- code, which is designed with 
tremendous mesh flexibility to handle a variety of grids 
including fully unstructured and non-orthogonal ones. 

2.1. Implementation within Cosmos+ + 

When working with more traditional spherical-polar 
meshes, the Cosmos-|— I- code actually evolves the 
MHD equations in a generalized coordinate system 
{xo, xi,X2, x^}, with the curvature implemented through 
metric terms. This is done even for the Newtonian for- 
mulation. The corresponding physical coordinates in 
general relativity are the Kerr-Schild polar coordinates 
{t,r,6,4}- For the cubed-sphere, instead, we construct 
the grid in physical space using the Kerr-Schild Cartesian 
coordinate system {t, x, y, z}. The two Kerr-Schild coor- 
dinate systems are related through the following trans- 
formations: 



x = r sm f cos cp — a sm u sm < 
y — r sin 9 sin 4 + d sin 6 cos < 
z = r cos , 



(1) 



{x^ + y'^ + - a^) + yj{x'^ + y2 + ^2 _ ^,2)2 _^ 402^2 
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2 , 2\ 1/2 

X +y^ 



cosO- 



cos ( 



ry — ax 
rx + ay 



v/(r2-Ha2)(x2-hy2) 

Ultimately the Cosmos-|--|- code just needs to know the 
coordinate locations of all the zone vertices. From those 
it is able to fully reconstruct all of the necessary zone 
properties such as volumes and face areas. We find it 
easiest for the cubed-sphere grid to start from the cubed- 
sphere coordinates {r, 41,42} of each vertex, then use 



(2) 
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the transformations given in Appendix [A| to convert to 
the Kerr-Schild polar coordinates {r, 6*7^, and finally 
use equation ([ij to obtain the correct Kerr-Schild Carte- 
sian coordinates {x, y, z} from the polar ones. For conve- 
nience we label the six blocks 0-5, with their orientations 
described in Appendix [X] Samples of blocks 0, 1, and 2 
are illustrated in Figure [IJ 

One consequence of using Kerr-Schild Cartesian coor- 
dinates is that, whereas in the spherical-polar case we 



Fragile et al. 


2007a|b|) 


ack hole musi 


J remain 



( Fragile fc Anninos 2005 
cubed-sphere case the b 
with the grid for a rotating black hole (a ^ 0). This is 
because in Kerr-Schild Cartesian coordinates, the z-axis 
is chosen to be the spin axis of the black hole and the 
event horizon is only symmetric about this axis. Thus, in 
order to get the inner boundary of the cubed-sphere grid 
to align with the black-hole event horizon, the black-hole 
spin axis must align with the grid z-axis. 

2.2. Boundary Conditions 

As we said, one of the difficulties with the cubed sphere 
is that the coordinates are not continuous across block 
boundaries. This requires some care when setting up 
communications between blocks. Even once the com- 
munication pattern between blocks is established, there 
are certain subtleties about the grid that must be dealt 
with. For instance, as we show in S|3] the gradient op- 
erator can only be made to converge properly (second 
order) if a set of ghost zones are constructed that are an 
extension of the coordinates on the current block. How- 
ever, such ghost zones then do not correspond directly to 
any of the zones on the neighboring block; instead they 
tend to straddle more than one zone, and a simple do- 
main exchange is not exactly valid. Fortunately, at any 
inter-block boundary it is only one of the (j) coordinates 
that is discontinuous; the radial coordinate and the other 
(f> coordina te are consistent acr oss any given inter-block 
boundary (Ronchi et al. 19961. Therefore, for a single 
layer grid with uniform zone spacing, the ghost zones of 
one block will never overlap more than two zones on the 
neighboring block. In such a case, we can get away with 
applying a boundary condition that simply fills the ghost 
zone with a field Fw that is a weighted average of the 
fields in the two real zones it overlaps, Fq and Fi. The 
weighted average we use is 

P Fo(L-|xo-x|) + fi(£~|xi-x|) 



where 



L = |xi - x| + |xo - x| 



(4) 



and x, xq, and xi are the coordinate centers of the ghost 
zone and the two real zones it overlaps on the neighbor- 
ing block, respectively. This weighting scheme is applied 
any time a normal domain exchange would be needed 
between neighboring processors, such as after fields are 
updated, but before any gradients are taken. Something 
slightly different must be done for advection as explained 
in the next section. 

2.3. Advection 

Because Cosmos++ was written using finite volume 
methods, and designed for arbitrary mesh topologies, few 



changes were needed to apply the code to the cubed- 
sphere mesh. The one thing (in addition to the ghost 
zone construction) that was modified, if only slightly, was 
the algorithm for advection. A number of different op- 
tions for advection are available in Cosmos+-f including 
upwind subzonal polyhedral reconstruction a nd global 
monotonic fl ux limiter methods, described in [Anninos] 



et al. (20051. These methods are designed to operate 
on multi-dimensional vector quantities (e.g., gradients) 
constructed from the convex attributes of arbitrary co- 
variant cell geometries and connectivities. However, for 
the cubed-sphere mesh we found that flux estimates per- 
formed with local one-dimensional limited projections (or 
differences) across individual cell faces are generally more 
robust than computing vector fluxes across the entire cell 
structure, even with appropriate multi-dimensional flux 
li miters. 

The me thod is only sl i ghtly modified from that pre- 



sented in Anninos et al. (20051, so we present only an 
abbreviated discussion. Fhe advection terms are solved 
for each evolved field quantity using an upwind time- 
explicit, first order forward Euler scheme with appropri- 
ately time-centered fluxes. Letting F represent any of the 
evolved flelds (or their consistent transport counterparts 
with F F/D, where D is the mass density), the dis- 
crete flnite- volume representation of the advection source 
term can be written 



faces 



(5) 



where Vz is the local donor cell volume of zone z, (Ai) f is 
the inward pointing area normal vector associated with 
face / of the donor cell, and (V) / is the face-centered 
velocity defi ned as a weighted aver age across neighbor- 



ing cells. In Anninos et al.| (|2 005| , the quantity (F*) 



, , , / 

represents piecewise linearly reconstructed zone-centered 
fields extrapolated to each cell face by a monotonic Tay- 
lor's series expansion, F* = F^ -I- {diF)^{r'^ — r\), pro- 
jected from the donor cell center to either the face 
center r* = r^- or the advection control volume center 
= r^. — (At/2)(y*) f, over a time-step interval Ai. The 

zone-centered limited gradient {diF)^ forces monotonic- 
ity in the extrapolated fields using polyhedral subzonal 
interpolations and control volume integrals to construct 
upwind, downwind and centered variations. The differ- 
ence here, for the cubed-sphere, is that the monotonic 
multi-dimensional gradient is replaced by a local one- 
dimensional calculation separately across each donor cell 
face, and along the direction of the cell face normal (per- 
pendicular to the cell face) using the generalized minmod 
limiter in the form 

1/a h\ 1/a c\ 1/6 

(6) 

where a ^ {I + A)VF£,, 6 = VFc, c = [l + A)VF[/, 
A is an order parameter between zero and unity spec- 
ifying the steepness of the applied limiter, and VFjy, 
VF^), and VFp are the upwind, downwind, and center- 
difference gradients, respectively. The upwind and down- 
wind gradients are defined as VF[/(£)) — k SF/Ss, where 
fc = ±1 depending on the upwind direction with respect 



VF, 



to the coordinate orientation, 6F 



is the 



Fig. 1. — Examples of blocks (right), 1 (center), and 2 (left) that might make up a cubed-sphere grid. Note that in this illustration we 
use a very low resolution for clarity. 



difference between donor and opposite cell field values, 
and Ss is the magnitude of the distance between donor 
and neighbor cell centers. The center-difference repre- 
sentation of the gradient is approximated as VFc — 
^j,^^g^(/c/2)((5F/(5s), where the sum is over opposite cell 
face pairs. A projected estimate for the advected fields 
contributing to the flux in equation ^ at each cell face 
is provided by the donor cell as F* = F^ -I- (5F = F^ — 
k \/F±xdr, where 6r = \xf-X:,-0.5At{V^Ai)A/{AiAj)\ 
is the distance to the advection control volume center 
along the direction aligned parallel to the cell face nor- 
mal vector A (between neighbor zone centers). 

For advection from one block to another, in order to 
conserve mass, energy, and momentum to round-off in- 
stead of truncation, it is important not to interpolate 
values between ghost zones as was done for the extrapo- 
lated field gradients in the previous section. Instead, for 
advection we use the ghost zones as "buckets" to cap- 
ture material advecting off of the host block. The mass, 
energy, and momentum collected in this bucket is then 
deposited into the corresponding real zone on the neigh- 
boring block that shares a face with the originating real 
zone as part of a final loop in the advection routine. This 
is appropriate since zones along inter-block boundaries 
share faces with only a single neighbor. 

3. GRADIENT TEST 

Because the cubed-sphere grid uses the same basic gra- 
dient operators that were already tested in Cosmos-!- + 
(Anninos et al. 20051, we fully expect the same second 
order convergence tor smooth fields, at least in the inte- 
rior zones. Nevertheless, it is worthwhile to conduct a 
simple gradient test for a variety of fields to verify sec- 
ond order convergence over the entire domain, including 
at the inter-block boundaries where we have introduced 
a new procedure for interpolation of fields beyond local 
grid domains. 

In our first attempt at implementing the cubed-sphere 
grid, we actually did not achieve uniform second or- 
der convergence. In that attempt, instead of construct- 
ing the ghost zones as extensions of each block as de- 
scribed in ^ we constructed the ghost zones to be ex- 
act replicas of the nearest zone on the neighboring block 



and to mimic the behavior of periodic boundaries on 
spherical-polar grids. However, this introduces a discon- 
tinuity into the gradient operator and actually prevents 
the convergence of gradients at the inter-block bound- 
aries. For interior zones not touching an inter-block 
boundary, we found the Ll-normalized error for the gra- 
dient of a simple scalar field to converge at second or- 
der as expected (the Ll-normalized error is defined as 
-^1 = Y^i,j,k\"-h3.k - A,],k\/{n^n.jnk), where at^.k and 
■^i,j.k E^re the numerical and exact solutions, respectively, 
in each zone and rii, nj, and rik are the number of zones 
in each of the three directions). However, for the interior 
zones touching the inter-block boundaries (not the ghost 
zones themselves, but the zones that touch them), the 
Ll-normalized error did not converge. To explain where 
this failure arises we first note that the gradient of a 
generic field F in Cosmos-|— I- is calculated as (akin to 
equation [s]) 

faces 

G.^d,F^--J2{F*A,)f, (7) 
^ / 

where the summation is performed over all cell faces. The 
problem arises in calculating F* , the face-centered value 
of the field; Cosmos-|— I- uses a simple average of the zone- 
centered values Fz in the two zones adjoining at face /. 
However, when the line connecting the two zone centers 
does not pass through the center of the zone face, as is 
the case for nearest neighbor cells across an inter-block 
boundary, this simple averaging does not give the correct 
face-centered value F* . In fact, it is relatively easy to 
show in this case that the absolute error {{aij^k ~ ^ij,fc|) 
in each zone along the inter-block boundary remains es- 
sentially constant, regardless of the resolution (it only 
depends weakly on the location of the zone along the 
boundary) , thus explaining the non-convergence in these 
zones. 

The ghost-zone construction described in fj2] on the 
other hand, which is the only one used for the remainder 
of this work, restores 2nd order convergence in all in- 
terior zones by giving a properly extrapolated value for 
F*. Here F* = 0.5{Fz -\- Fyy) is a simple average of the 
zone-centered value Fz in the interior zone and the ghost- 
zone weighted average Fw from equation ([3|. We have 
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confirmed that all interior zones (including those touch- 
ing the inter-block boundaries) give errors at the level of 
round-off for flat fields and second-order convergence for 
all linear and higher-order fields. 

4. TILTED ACCRETION DISKS 

Having demonstrated that our implementation of the 
cubed-sphere grid preserves the correct convergence or- 
der for our code, we can confidently move on to testing 
our primary application of interest: black-hole accretion 
disks. We begin with a review of how the simulations 
are init ialized and then consider two sets of test cases: 
In §4.1| we study disks of diffe ring alignments relative to a 
Schwarzschild black hole; in [4.2 we compare tilted disk 
simulations around a Kerr black hole, one carried out 
on a spherical-polar mesh and the others on the cubed- 
sphere. 

Most of the accretion disk simulations presented in this 
work using the cubed-sphere grid are carried out at a res- 
olution of 128 X 64 X 64 X 6, where there are 128 radial 
shells and each of the blocks are resolved with 64 x 64 an- 
gular zones. Along its symmetry planes, such a grid looks 
like a spherical polar grid of resolution 128 x 128 x 256. 
However, the more uniform distribution of zones in the 
cubed-sphere grid means we are able to achieve such reso- 
lution with a smaller number of zones overall (by a factor 
of 3/4). Also, because of the more uniform zone sizing, 
we are able to run with a Courant time step that is al- 
most 30 times larger than could be used with a spherical- 
polar grid of that resolution, which means the required 
CPU cycle-count is smaller by the same factor. An image 
of the actual grid used in these simulations is shown in 
the left panel of Figure [2j This can be compare d to the 
spherical-polar grid used in [Fragile et al. (2007b I, includ- 
ing the underresolved polar regions, which is shown in 
the right panel of Figure [2j The timestep for the cubed- 
sphere grid is even 25% larger than for that special grid, 
where the pole was underresolved precisely to keep the 
timestep reasonable. 

The inner and outer radial boundaries are set at 
0.98rBH and 120rG, respectively, where tbh is the ra- 
dius of the black- hole horizon and re = GM-qh/c^ is the 
gravitational radius. Note that, because we use the Kerr- 
Schild form of the Kerr metric, we are able to place the 
inner radial boundary inside the black-hole horizon. In 
principle, this should keep the inner boundary causally 
disconnected from the flow, although numerically there 
is still some communication. At both the inner and outer 
radial boundaries we apply "outflow" conditions: Fluid 
variables are set the same in the external boundary zone 
as in the neighboring internal zone, except for velocity, 
which is chosen to satisfy 



mt 
* int 



y points off the grid , 
y points onto the grid . 



(8) 



For the initial conditions of the simulations we start 
from the commonly used analytic solution for a hydro- 
static fluid torus orbiting the black hole. In this case, 
we choose the torus parameters to be the black-hole spin 
(o/Mbh), the inner radius of the torus (rjn = 15rc), the 
radius of the pressure maximum of the torus (reenter — 
25rG), and the power-law exponent (q = 1.68) used in 
defining the specific angular momentum distribution, 

e = -u^/ut - fcA2-9 , (9) 



where = g^^u'^ , g^^ is the 4-metric, and is the fiuid 
4- velocity. We then follow the procedure in |Chakrabarti| 
(fl985| to solve for the initial state of the torus. Knowl- 
edge of reenter Icads directly to a determination of Reenter 
by setting it equal to the geodesic value at that radius. 
The numerical value of k comes directly from the choice 
of q and the determination of Acentcr, where 



1 

A2 



5t0 + ^9tt 



Finally, having chosen r^. 



,2 • (10) 

we can obtain Um = ut{rin), 
the surface binding energy of the torus, from u^^ — g** — 

2V'^-|-^V^- 

The solution of the torus variables c an now be speci- 
fied. The internal energy of the torus is ( De Villiers et al 
2003| 



ut{r,e)f{e{r,e)) 



(11) 



where £in — ({rm) is the specific angular momentum of 
the fiuid at the surface and 



1 - fc2/»^" 



1/ct 



(12) 



where n = 2 — q and a — (2n — 2)/n. Assuming an isen- 
tropic equation of state /or the initialization only, the gas 
pressure and density must be related by the expression 
P = pe(r — 1) — Kp^ , and so the density is given by 

p = [e(r - 1)/k]^/*^"^^ We take F = 5/3 and n = 0.01 
(arbitrary units). Finally, the angular velocity of the 
fiuid is specified by 



y4, ^ 9t,p + igt: 



94><t> + ^9t<t> 



(13) 



The torus is then seeded with weak poloidal mag- 
netic field loops with non-zero spatial components = 



d{,A^ and ~ drA^, where 



A, 



b{P - Pcut) for p > Pent 
for p < Pent 



(14) 



The parameter p^ut — 0-5*Pniax,o is used to keep the field 
a suitable distance inside the surface of the torus initially, 
where Pmax.o is the initial density maxi mum within the 
torus. Using the constant b in equation (14), the field is 
normalized such that initially /3mag = P/Tb ^ /5mag,o = 
10 throughout the torus, where Pb is the magnetic pres- 
sure. The magnetic field is added in order to seed the 
magn eto-rotational instability (MRI; [Balbus fc Hawley 
1991 1, which is now commonly believed to be the source 



of angular mom entum tran sport within black-hole accre- 
tion disks (Balbus & Hawley 1998). 

As mentioned in ij2] our implementation of the cubed- 
sphere requires that the black-hole be aligned with the 
grid. Therefore, unlike our previous work where we tilted 
the black hole, if we want a tilted configuration now we 
must tilt the disk. By itself, tilting the disk is a rather 
trivial operation, simply requiring the following coordi- 
nate transformation be applied prior to constructing the 
torus: 



x' = x cos f3o — z sin /?o 

y'^y 

z' = x sin Po + z cos /3o , 



(15) 
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Fig. 2. — {'e/i) Plot of the cubed-sphere gr id geometry used for the disk simulations presented in this work, (right) Plot of the 
spherical-polar grid used in [Fragile et al.| | |2007b| l, including an underresolved polar region. 



where f3o is the initial tilt of the disk. However, as we de- 
scribe in the next section, we were surprised to discover 
that our tilted disks evolved differently than our un- 
tilted disk, at least at early times, even for non-rotating, 
Schwarzschild black holes, for which a tilt should have 
no physical meaning or significance. 

4.1. Schwarzschild Black Hole 

Here we compare simulations of black-hole accre- 
tion disks carried out for a Schwarzschild black hole 
(o/Mbh = 0). Table [l] summarizes the parameters for 
these runs. In our naming convention, the first number 
indicates the dimensionless spin of the black hole (a/M) 
without the decimal; the second number, if present, gives 
the tilt angle of the disk in degrees; the final letter is used 
to distinguish what resolution the simulation is carried 
out at, "H" being our high resolution (128 x 64 x 64 x 6) 
and "L" being low (64 x 32 x 32 x 6). The two main 
simulations, OH and 015H, begin from identical initial 
conditions except for the tilt of the disk with respect 
to the grid, which are and 15° respectively. We also 
include results of a simulation that uses the spherical- 
polar grid from Fragile et al. (2007b I; this simulation is 
denoted by the suffix "SP" arid h as an equivalent peak 
resolution of 128^. We showed in Fragile et al. (2007b I 
that this was roughly the minimum resolution needed to 
get a relatively well converged result for this type of prob- 
lem. Therefore, in the current work, we do not expect 
our low-resolution simulation (OL) to be converged; they 
are instead included for the purpose of estimating the 
rate of convergence when using the cubed-sphere grid. 

Our first concern with the cubed-sphere grid is that 
the angular-momentum conservation may not be suffi- 
cient for the purpose of following the long-term evolu- 
tion of an accretion disk, particularly as the flow crosses 



the coordinate discontinuities at block boundaries. At a 
minimum we want to quantify our angular-momentum- 
conservation error, which we do graphically in Figure [3] 
where we plot the total angular momentum in each sim- 
ulation as a function of time for runs OL, OH, 015H, and 
OSP. The total angular momentum is defined as 

(16) 

Jv 

where = {ph + 2Pb)u°u^ - {B^B^)/{4tt), g is the 
determinant of the 4-metric, B^ is the magnetic field 4- 
vector, and 



h = 1 



P 
P 



(17) 



is the relativistic enthalpy. We only plot part of the first 
orbital period (t < 0.8torb) of data because after this 
time significant amounts of angular momentum begin to 
advect into the black hole and leave the outer boundary 
of the grid in jets and winds, so that it is much more 
difficult to track the global conservation. Ideally all the 
lines in Figure [3] would be perfectly horizontal, indicat- 
ing exact angular momentum conservation, but we do 
not really expect this (slight imbalances in the momen- 
tum "source" terms and imperfect boundary conditions, 
for instance, can prevent exact conservation). The an- 
gular momentum conservation in our "low" resolution 
simulation OL is 0.18% (extrapolated to a full orbital pe- 
riod); this drops down to 0.048% per orbital period at our 
normal resolution, about the level of convergence (second 
order) we expect. Furthermore, it appears that this error 
is not strongly dependent on the orientation of the disk 
with respect to the grid, based on a comparison of simu- 
lations OH and 015H. Simulation OSP is included to give 
some indication of our typical angular momentum conser- 
vation error on the multi-resolution-layer spherical-polar 
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TABLE 1 

SCHWARZSCHILD SIMULATION PARAMETERS 



Simulation 


a/M 


Tilt 


Resolution^ 


Start^ 


Endb 






Angle 




Time 


Time 




0L<1 








64 X 32 X 32 X 6 





4 


-0.0129 










128 X 64 X 64 X 6 





4 


-0.0090 


015Hd 





15° 


128 X 64 X 64 X 6 





4 


-0.0085 


OSP<= 








128^ 





4 


-0.0141 



In the case of the spherieal-polar grid this represents the equivalent peak resolution of an unrefined grid. 
In units of toj-b — 785Gi\//c^, the geodesic orbital period at the initial pressure maximum reenter- 
Calculated from the slopes of M vs. t over the interval 3 < t/to^h ^ 4. 
Cubed-sphcre grid. 

Multi-resolution-layer spherical-polar grid. 




Fig. 3. — Plot of the total angular momentum as a function of 
time for simulations OH {solid), 015H (dashed), OL (dot-dashed), 
and OSP (dotted). All plots have been normalized to their initial 
angular mo menta. Simulation OS P uses the spherical-polar grid 
described in [Fragile et al.|p007b[ |. 

grid used in our previous work. The error in this case 
is 0.019% per orbital period, somewhat better but still 
comparable to simulation OH, suggesting we suffer only a 
small degradation in angular momentum conservation in 
going from our spherical-polar grid to the cubed-sphere 
grid. 

Because we are simulating a non-rotating black hole in 
this section, any tilt we assign the disk has no physical 
meaning; it can only be defined relative to the grid. We 
would expect, therefore, that this tilt would not have any 
physical efi^ect on the evolution. Interestingly, that is not 
what we find at early times. The difference is perhaps 
shown most graphically in Figure [4] which shows the gas 
density of the disk for simulations OH and 015H along one 
azimuthal slice after one orbital period (torb) at the initial 
pressure maximum (reenter)- In simulation OH (left panel 
of Fig. |4| , the disk has spread radially to such an extent 
that it reaches all the way to the event horizon of the 
black hole (inner boundary of the computational grid) . 
In simulation 015H, on the other hand {right panel), the 
disk has hardly spread radially at all, having started at 
Tin = 15rG and only penetrated to 12rQ. 

The primary mechanism responsible for the radial 
spreading of the disk over the first orbital period is not 
solely the MRI, but also the differential winding of the 
initial radial component of the poloidal field loops, the 
so-called il-dynamo. The amplified toroidal and radial 
field components allow for efficient angular momentum 
transport essentially from the beginning of the simula- 
tion. This is, of course, peculiar to an initial field con- 



figuration such as ours which includes a radial field com- 
ponent. If, instead, we started from a purely toroidal 
field, differential winding would not play a role initially 
and angular momentum transport would have to await 
a more complete development of the MRI, which occurs 
on roughly an orbital timescale. 

Something in simulation 015H appears to be shorting 
out the shear amplification of the field as compared to 
simulation OH. Growth of the MRI also appears to be 
delayed, as evidenced by the less turbulent appearance 
of simulation 015H in Figure |4] This may be related to 
the lack of an fi-dynamo since the MRI has less field 
to gr ow on whenever this is inactive (Hawley & KrolikI 
20021. Furthermore, we can see for certain in J^'igure [s] 
that the total magnetic energy is growing more slowly in 
simulation 015H than in OH (and OSP). Here we define 
the magnetic and kinetic energies as 



(ff°° + 2uV)PB 



47r 



(18) 



and Dh{vP — 1), respectively, where D = W pis the gen- 
eralized fluid density with boost W = ^J—gvP. Both 
energies are summed over the entire simulation domain. 
All three simulations show a very rapid initial growth 
of the magnetic energy due to the combination of shear 
amplification and the MRI. They also show a gradual 
increase in kinetic energy over the first orbit as gravita- 
tional potential energy is converted into kinetic. After 
approximately Itorb the growth of the magnetic fields 
saturates. At about the same time in simulations OH 
and OSP, kinetic energy begins accreting into the black 
hole in significant amounts, accounting for the sudden 
change in slope. This happens about an orbit and a half 
later in simulation 015H. 

The culprit for the retarded field growth in simulation 
015H appears to be the numerical treatment of the cur- 
rent sheet that forms in the midplane of the disk as a 
result of the differential winding. For an untilted simu- 
lation, such as OH, this current sheet resides almost ex- 
actly along an interfacial boundary, right along one of the 
symmetry planes of the grid (see left panel of Figure |6| . 
Furthermore, because this is a nearly perfect symmetry 
plane for the flow, there is very little advection of fluid 
across this boundary, and so the current sheet remains 
relatively stationary. In effect, the current sheet remains 
unresolved, because it spans less than a full zone's width 
in the vertical direction. Notice how narrow the cur- 
rent sheet is in the left panel of Figure |6] This is not 
the case for a tilted-disk simulation. By necessity, the 
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Fig. 4. — Plot of logarithm of density (normalized to po,max) along an azimuthal slice at i/i = at t = Itorb for simulations OH {left 
panel) and 015H (right panel). Because the black hole is not rotating in this simulation, the tilt should have no physical effect and we 
would expect the two simulations to evolve nearly identically. The observed differences are due to the numerical treatment of the current 
sheet that forms in the midplane of the disk, as described in the text. 

sheet forms if one starts from a purely toroidal field. Fur- 
thermore, as the disk becomes more turbulent with the 
action of the MRI, we find that the discrepancies between 
the tilted and untilted simulations are dramatically re- 
duced to the point that, at late times, they are nearly in- 
distinguishable. For instance, in Figure [7] we show plots 
equivalent to Figure [4] except a,t t = Atoih as opposed to 
l^orb, which show the two disks to be nearly identical. 
The late-time mass accretion rates are also quite similar 
(see Table [T|. 

For a more rigorous comparison, in Figure |8] we 
present time- and shell-averaged values of density (p), 
gas pressure (P), dimensionless stress (a), plasma mag- 
netization parameter (/3mag), specific angular momentum 
{£), and radial inflow velocity {V^) as functions of radius 
for simulations OH, 015H, and OSP at late time, where 
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Fig. 5. — Plot of the total magnetic (fop) and kinetic (bot- 
tom) energies as functions of time for simulations OH (solid), 015H 
(dashed), and OSP (dotted). All plots have been normalized by the 
initial kinetic energy of simulation OH. 



disk midplane in a tilted disk is no longer aligned with 
any symmetry plane of the grid (fright panel of Figure 
[g]). This means that the disk midplane, and more impor- 
tantly the midplane current sheet, passes through the 
interiors of some zones rather than always along their 
boundary. Numerically this is a critical distinction. For 
a zone-centered code such as CosmosH — h whenever a cur- 
rent sheet is aligned along an interfacial boundary that 
experiences no advection, as is approximately the case 
in our untilted simulations (OH and OSP), there can be 
no numerical reconnection and magnetic fields are pre- 
served. If, on the other hand, the current sheet passes 
through a zone center, as it does in our tilted simulation 
(015H), numerical reconnection is greatly enhanced. The 
effect is to drain energy from the magnetic field. In the 
present work, which uses the internal energy conserving 
form of Cosmos -) — h this energy is simp ly lost from the 
simulations (see Fragile & Meier (20081) for a discussion 
of the implications of the different forms of energy con- 
servation in numerical simulations of black-hole accretion 
disks) . 

This is a somewhat worrisome discovery; however, we 
emphasize that it is restricted to the particular field ge- 
ometry we start with, as no strong midplane current 



\vJ'u^\\B\\'^ - B''B'^\ 



(19) 



and V — {pV^)/{p). Angle brackets indicate that a 
radial shell-average has been taken, where 

{Q){r,t) - ^ r r Q^dedcj, , (20) 



1 

A 

and A — J^^ V~~9 dcj) is the surface area of a given 
radial shell. The time-averaging is done over the interval 
3 < i/torb < 4. The shell-averages for P, a, /3mag, ^, and 
V are mass- weighted. Measurements of p, P, and £ show 
very good agreement between all three simulations, with 
errors everywhere < 20% and generally much less. The 
discrepancies in a, /3mag) and V are similarly small for 
simulations OH and OSP, but considerably larger for sim- 
ulation 015H. This is not unexpected as these quantities 
depend sensitively on the distribution of magnetic field, 
meaning they are more affected by the delayed growth of 
the MRI. 

4.2. Kerr Black Hole 

Having shown that the late-time evolution of simulated 
black- hole accretion disks on our cubed-sphcrc grid is rel- 
atively independent of the orientation of the disk with 
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Fig. 6. — Pseudocolor plot representing the value of (in code units) along an azimuthal slice at </> = at t = Itorb for simulations OH 
{left panel) and 015H {right panel). The midplane current sheet (represented by the line where the color changes from red to blue) remains 
essentially sub-zonal in simulation OH, whereas it is spread across approximately 3 zones in simulation 015H (greenish-yellow zones between 
red and blue). 
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Fig. 7. — Same as Fig. |4] except at time t = it^^^, instead of Itorb- 
{right panel) simulations are greatly diminished. 

respect to the grid by analyzing a few Schwarzschild test 
cases, we can now evaluate the treatment of tilted accre- 
tion disks around modestly rotating (o/Mbh = 0.5) Kerr 
black holes. Here our test simulations (515L and 515H), 
which use the new cubed-sphere grid, are compared to 
a reference simulation (515SP), which uses the multi- 
resolution-layer spherical- 
( 2007b I (shown in Figure [2] right panel) 
angle 



olar grid from Fragile et al. 

All simulations 

have an initial tilt angle ~3o — 15°. Again we do not 
expect the low-resolution simulation (515L) to be con- 
verged; instead it is included to provide some indication 
of the rate of convergence. The parameters for each run 
are described in Table |2] 
First we show in Figure [9] that the general disk prop- 



0.0085 
Max: 0.8096 
Min: 8.972e-10 



Here the discrepancies between the untilted {left panel) and tilted 

erties of simulations 515H and 515SP are quite similar. 
Again, the largest discrepancies are in the dimension- 
less stress a and the plasma magnetization parameter 
Anag = P/Pb- This is not surprising since both of these 
properties have been show n in previous studies to be 
quite sensitive to res olution (Hawley et al. 1996 Fromang 



&: Papaloizou 20071, and although the total number of 
zones in these two simulations is comparable, the distri- 
bution of those zones is considerably different. The level 
of agreement in the other parameters is really quite re- 
markable given the very different structures of the grids. 
Of course, this was exactly what we were hoping to see. 

Now, because the black-hole is rotating, the tilt of the 
disk has some physical meaning and consequently causes 
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Fig. 8. — Main disk properties plotted as a function of radius for simulations OH, 015H, and OSP. The data have been time-averaged over 
the final orbital period of each simulation (3 < t/torb 5; 4). P, a, /3, £, and V" are mass-weighted averages of the pressure, dimensionless 
stress, plasma equipartition parameter, specific angular momentum, and radial inflow velocity, respectively. 
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TABLE 2 
Kerr Simulation Parameters 



Simulation 


a/M 


Tilt 


Resolution^ 


End*^ 






Angle 




Time 




515Ld 


0.5 


15° 


64 X 32 X 32 X 6 


10 


-0.0032 


515H<i 


0.5 


15° 


128 X 64 X 64 X 6 


10 


-0.0114 


SISSP'^ 


0.5 


15° 


128^ 


10 


-0.0122 



^ In the case of the spherical-polar grid this represents the equivalent 
peak resolution of an unrefined grid. 

^ In units of ^orb ^ 789GM/c^, the geodesic orbital period at the 
initial pressure maximum reenter. 

Calculated from the slopes of M vs. t over the interval 3 < t/torb ^ 4. 
^ Cubed-sphere grid. 

^ Multi-resolution-layer spherical-polar grid. 
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Fig. 10. — Plot of the tilt /3 as a function of radius through the 
disk for simulations 515L, 515H, and 515SP. The data for this plot 
has been time averaged over the final two orbital periods of each 
simulation (8 < t/torb 10). The initial tilt was /3o = 15°. 



as 
Blaes' 



changes in i ts evolution re l ative t o an untilted disk, 
describ ed in 'Fragile et al. ( 2007b I and Fragile 
(|2008|. For instance, although tne disk begins with a 
uniform tilt of /?o = 15° , we expect a warp caused by the 
gravitomagnetic torque of the black hole to propagate 
through the disk as a bending wave. This will cause the 



tilt (3 to become an oscillating furiction of radius ( Ivanov 



fc Illarionov l 11997 
we plot [5{r 



Lubow et al. 20021. In Figure TU 



time averaged over the interval St, 



orb 



< 



t < lOiorb, for simulations 515L, 515H, and 515SP. As 



in prev ious work ( Fragile fc Anninos||2005 Fragile et al. 
2007b| 



we recover the tilt using the dehnition 

JbH • JDisk(?') 



arccos 



'BHlHDisk 



{r)\ 



(21) 



where Jbh is the angular momentum vector of the black 
hole and JDisk('') is the angular momentum vector of 
each radial shell of the simulation domain (dominated 
by the disk). Again simulations 515H and 515SP pro- 
duce remarkably similar results, with discrepancies no 
larger than 10% and generally much smaller. The dis- 
crepancies likely have their root in the small differences 
in conditions at the inner edge of the disk (see Figure ^ 
where the bending waves are launched. The 515L simula- 
tion exhibits considerably larger discrepancies over most 
of the disk as expected. 

Along with warping the disk, the gravitomagnetic 
torque of the black hole also causes it to precess, partic- 
ularly in disks such as the ones in our simulations where 




Fig. 11. — Plot of the precession (twist) 7 as a function of time 
for simulations 515L, 515H, and 515SP. The slope of this plot can 
be used to estimate the precession period of the disk as a whole, 
which is 0.7(M/Mq) s. 

the fast sound-crossing time causes the disk material to 
be tightly coupled and respond globally to the torque of 
the black hole. Global precession of this nature h as been 
noted b efore in low Mach number hydrodynamic (Fragile] 
fc Anninos 2005) and MHD ( Fragile ct al. 2007b| ) disks. 
We track tke overall precession (twist), defined as 



7 ^ 



'BH 



X J 



Disk 



'BH 



X J 



Disk I 



y 



(22) 



where Joisk is the total angular momentum vector of the 
disk and y is the unit vector that points along the ini- 
tial line-of-nodes between the black-hole symmetry plane 
and disk midplane. In order to capture twists larger than 
180°, we also track the projection of Jbh x Joisk onto x, 
allowing us to break the degeneracy in arccos. By plot- 
ting the cumulative precession as a function of time as we 
have done in Figure [TT] we make it easy to calculate the 
precession period of the disk - in this case O.7(M/M0) s, 
which agree s nicely with our pred ictions for a black-hole 
of this spin (IFragile et al.||2007bl. 



5. CONCLUSION 

In this paper we have presented our implementation 
of the cubed-sphere grid within Cosmos-|— 1-. The cubed- 
sphere grid has at least three significant advantages over 
more-traditional grid options: 1) it has topological prop- 
erties similar to a Cartesian grid, but generally conserves 
angular momentum much better (and nearly as well as a 
spherical-polar mesh); 2) it can run at a larger Courant- 
limited timestep than a spherical-polar mesh at compa- 
rable resolution (almost a factor of 30 at the resolution 
used in this work); and 3) it distributes zones more evenly 
than a spherical-polar mesh, which is desirable for prob- 
lems where the symmetry is imperfect, such as in tilted 
accretion disks around rotating black holes, a problem of 
particular interest to us. 

In Section[2]and Appendix|X]we gave detailed prescrip- 
tions for the construction of the cubed-sphere grid and 
shared a few "lessons learned" in regards to extrapolat- 
ing fields at inter-block boundaries and applying limiters 
to the field gradients in advection. After implementing 
these lessons ourselves, we found we could recover second 
order convergence over the entire grid, including along 
the inter-block boundaries. 

To specifically demonstrate that the cubed-sphere grid 
is a viable option for the black-hole accretion disk work 
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we have in mind, we have carried out a series of such sim- 
ulations on our new cubed-sphere grid, using results from 
our spherical-polar grid as a reference standard. From 
these tests we conclude that: 

• The cubed-sphere grid conserves angular momen- 
tum nearly as well a spherical-polar grid at compa- 
rable resolution. 

• The angular momentum conservation error on the 
cubed-sphere grid is only weakly dependent on the 
tilt of the disk. 

• Results on the cubed-sphere grid converge to the 
same solutions obtained on a spherical-polar grid 
when the two grids approach comparable resolu- 
tions. This is true for both untilted and tilted disks. 

• Important disk properties such as density, pres- 
sure, specific angular momentum, inflow velocity, 
tilt and twist agree to better than 10-20% for sim- 
ulations carried out on cubed-sphere and spherical- 
polar grids with roughly (2-3) x 10^ zones. 

During our testing, we made one surprise discovery 
- that the early-time evolution was considerably differ- 
ent between our untilted and tilted simulations on the 
cubed-sphere grid. We found this to be true even for 
non-rotating Schwarzschild black holes, for which a tilt 
should have no physical meaning or significance. This is 
something we had not seen on the spherical-polar grid, 
but there we had tilted the black-hole, not the disk as we 
do now. We did not anticipate how important this dif- 
ference would be for the initial growth and development 
of the fi-dynamo and MRI. 

We attribute the disparate early-time behavior to the 
differing ways in which the strong initial current sheet in 
our disk is handled numerically when it is tilted with re- 
spect to the grid. This is another reminder of the impor- 
tant role that numerical reconnection plays in the evo- 
lution of numerically simulated magnetized flows even 
though this topic is perhaps not given enough empha- 
sis in the literature. The appearance of current sheets 
is virtually unavoidable in strongly sheared MHD flows 
such as accretion disks. One possible technique for treat- 
ing the current sheets more consistently throughout the 
simulation may be to use an artiflcial resistivity. This 
would ensure that the current sheets are always resolved 
in a similar fashion regardless of their orientation with re- 
spect to the grid. However, this technique has o nly been 
implemented very recently in r elativistic MHD (K omis- 
[sarov 2007) . An alternative, although only partial, so- 
lution nught be to use a total-energy conserving scheme 



instead of the internal-energy conserving one used here. 
This, at least, guarantees that the energetics of the flow 
remain consistent by recapturing in the form of thermal 
energy any energy lost through magnetic reconnection. 
When coupled with a radiative cooling package, this can 
give a muc h more physical descri ption of the evolution 
of the flow ( [Fragile fc Meier||2008[ ). 

Although the numerical treatments of current sheets 
and reconnection are important to understand and ap- 
preciate, it is equally important in the context of this 
paper to point out that we demonstrated by numerical 
example that the long-term evolution of our disks is rel- 
atively unaffected by whether or not they are tilted with 
respect to the grid. As expected, only when the tilt is 
relative to a rotating black hole are there long-term im- 
plications within the disk. 

We are not surprised to find significant discrepancies 
between our "low" and "high" resolution simulations, as 
previous experience had shown us that 128'^ was roughly 
the minimum resolution necessary to follow the evolution 
of black- hole accretion disks in global general-relativistic 
MHD simulations such as these. Below that resolution 
the characteristic MRI wavelength (Amri = 2ttva/^, 
where va is the Alfven speed) is not covered by a suf- 
ficient number of zones over much of the disk volume. 
This has nothing to do with the cubed-sphere grid it- 
self, but is rather a universal constraint for these types 
of problems. 

Overall we consider our experimentation with the 
cubed-sphere grid to be a success. In future work we 
will present further analysis of tilted disks (and their as- 
sociated jets) evolved using this new grid option. 
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APPENDIX 

CUBED-SPHERE TRANSFORMATIONS 

Included in this appendix are the transformations necessary to go from the cubed-sphere coordinates {r, 0i , 
the corresponding spherical-polar ones {r, 9, ((>} on each block. 

• Block - (centered about the -l-x-axis; 7r/4 < 0i < 37r/4; — 7r/4 < 02 < 7r/4) 

cos (pi cos (1)2 



i} to 



cos 9 : 
sin 9 = 



\/l - (cos 01 sm4>2)^ 



sin = sin (j)2 

cos (j) = COS ^2 



Block 1 - (centered about the +y-axis; 7r/4 < (pi < 37r/4; 7r/4 < 4'2 < 37r/4) 

cos ^1 sin (p2 
cos C = . = 

— (cOS(?iiCOS(^2)^ 

sin 6 = \/l — cos^ ^ 

sin = sin (/)2 

cos (j) — cos (/)2 

Block 2 - (centered about the — x-axis; 7r/4 <4>\< 37r/4; 37r/4 < (j>2 < 57r/4) 

^ — cos COS (}!)2 

COS f) = ^= 

\/l - (cos 01 sin(/)2)2 

sin = \/l — cos^ 9 
sin (/> = sin <^2 
cos = cos 02 

Block 3 - (centered about the -?/-axis; 7r/4 < <pi < 37r/4; 57r/4 < 02 < 77r/4) 

^ - cos 01 sin 02 
cos = — 

— (cos 01 cos 02)^ 

sin 9 = — cos^ 9 

sin = sin 02 

cos = cos 02 

Block 4 - (centered about the +z-axis — 7r/4 < 0i < 7r/4; — 7r/4 < 02 < 7r/4) 

cos 01 cos 02 
cos w = = 
a/1 — (sin 01 sin 02)^ 

sin 9 = \/l — cos^ 9 

. , cos 01 sin 02 

sm = — 

sin^y'l — (sin 01 sin 02)^ 

sin 01 cos 02 
cos — — ^ ^ = . 

sin Oy'l — (sin 0i sin 02)^ 

Block 5 - (centered about the — ^-axis — 7r/4 < 0i < 7r/4; — 7r/4 < 02 < 7r/4) 

— cos 01 cos 02 

COS 9 = = 

i/l - (sin 01 sin 02)2 

sin 9 = \/l — cos^ 9 

. , cos 01 sin 02 

sm = — 

sin^y'l — (sin 0i sin 02)^ 

— sin 01 cos 02 

cos = = . 

sin6'\/l — fsintii smSo)"^ 
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